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Abstract 

In this paper we consider the stability and convergence of numerical 
discretizations of the Black-Scholes partial differential equation (PDE) 
when complemented with the popular linear boundary condition. This 
condition states that the second derivative of the option value vanishes 
when the underlying asset price gets large and is often applied in the ac- 
tual numerical solution of PDEs in finance. To our knowledge, the only 
theoretical stability result in the literature up to now pertinent to the lin- 
ear boundary condition has been obtained by Windcliff, Forsyth & Vetzal 
|14| who showed that for a common discretization a necessary eigenvalue 
condition for stability holds. In this paper, we shall present sufficient con- 
ditions for stability and convergence when the linear boundary condition 
is employed. We deal with finite difference discretizations in the spatial 
(asset) variable and a subsequent implicit discretization in time. As a 
main result we prove that even though the maximum norm of e tM (t > 0) 
can grow with the dimension of the semidiscrete matrix M, this generally 
does not impair the convergence behavior of the numerical discretizations. 
Our theoretical results are illustrated by ample numerical experiments. 
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1 Introduction 



A popular assumption in the valuation of financial options via the numerical 
solution of partial differential equations (PDEs) is the so-called linear boundary 
condition, see, for example, [TJ [TUJ [T31 [H] . The linear boundary condition states 
that the second derivative of the option value with respect to the underlying 
asset price vanishes if the asset price gets large. This condition represents a 
linear behavior of the option value for large asset prices, which can be seen to 
hold for a wide variety of financial options. In spite of its broad use in prac- 
tice, only few rigorous results have been derived in the literature up to now on 
the stability and convergence of numerical discretizations if the linear boundary 
condition is applied. As it turns out, in the finite difference (FD) approach a 
natural treatment of the linear boundary condition leads to a downwind dis- 
cretization of the advection term at the relevant grid point; the details of which 
are given below in this section. Consequently, in the actual numerical solution 
one might expect instability, or at least an adverse effect on the convergence 
behavior. It appears, however, that this is not observed in practice. To our 
knowledge, the only theoretical stability analysis in the literature up to now 
pertinent to the linear boundary condition has been performed by Windcliff, 
Forsyth & Vetzal [14]. These authors proved that for a common discretization 
of the Black-Scholes PDE a necessary eigenvalue condition for stability holds. 
The objective of the present paper is to arrive at useful sufficient conditions for 
stability and convergence of discretizations when the linear boundary condition 
is employed. As far as we are aware, such conditions are lacking in the current 
literature, but they are clearly of much interest. 
Consider the Black-Scholes PDE 

du d^u du 

— ( Sl t) = \a 2 s 2 -^ I {s 1 t) + rs — {s,t)~ru{ S ,t) (s > 0, < t < T), (1.1) 

where r > and a > are given real constants that denote the risk-neutral 
interest rate and the volatility, respectively, and T > is the given maturity 
time of the option. The exact solution u(s,t) represents the fair value of an 
option if the underlying asset price equals s at time T — t. 

For the numerical solution, one restricts in practice the s-domain to a boun- 
ded set [0,5] with fixed 5 > chosen sufficiently large. The PDE (JUJ is 
complemented with initial and boundary conditions. In this paper, we consider 
at s = 5 the linear boundary condition 

u ss (5,i) = (0<t<T). (1.2) 

At the lower boundary s — a standard Dirichlct condition is taken, which 
depends on the particular option. The initial condition is given by the payoff of 
the option. 

FD discretization of the initial-boundary value problem for (jl.ip on a general 
(non-uniform) grid = sq < s\ < S2 < ■ ■ ■ < s m +i < «m+2 = 5, with mesh 
widths hj — Sj — Sj_i, leads to an initial value problem for a system of ordinary 
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differential equations (ODEs), 



U'(i) = MU(t) + b(t) (0<t< T), U(0) = U ■ (1.3) 

Here M denotes a given real (m + 2) x (m + 2)-matrix and Uo and b(t), for 
< t < T, are given real (m + 2)-vectors. The vector Uo is directly given by 
the payoff function and b depends on the Dirichlct condition at s = 0. In this 
paper we shall deal with matrices M of the form 



M = 



A 


B ' 




C 



( OL\ 



V 



7i 
o 2 



72 



A, 



A, 



rS 

h m+2 
rS 



n+2 



(1.4) 



where a?, /3j, 7j denote given real numbers. The 2 x 2-matrix C represents a 
natural discretization of the linear boundary condition (|1.2[) . It is determined 
by the following approximations at the grid points s m +i and s m +2 = S: 



Us(s m +l,t) 
U s (s m +2,t) 



u(s m+2 ,t) - u(s m+ i,t) 

h m +2 

u(s m+2 ,t) - u(s m+ i,t) 
h m +2 



u ss (s m+ i,t) w 0, (1.5a) 
u ss (s m+2l t) = 0. (1.5b) 



As r in (|1.1[) is positive, the approximation of u s forms an upwind scheme at 
Sm+iy but the same approximation constitutes a downwind scheme at s m +2- The 
latter approximation can be regarded as obtained from the second-order central 
scheme for advection at s m+ 2 with virtual point s m+3 := s m+2 + h. m+2 and then 
replacing u(s m+ 3,t) by 2w(s m+2 ,i) — u(s m+ i,t) in view of the linear boundary 
condition. The discretization (jl.5b ) at s m +2 is identical to the one considered 
in |14j . The discretization (jl.5b ) at s m +i, on the other hand, appears to be 
new. In particular, we approximate u ss by zero at this point instead of using 
the standard second-order central scheme for diffusion. The choice ()1.5j) yields 
a partial decoupling between the FD solution at the grid points s m +i, s m +2 and 
that at si, S2, . . . , s m . Concerning the discretization on [si,s m ] we make no 
assumptions yet, except that at each relevant grid point Sj the stencil belongs 
to {sj-i, sj, Sj+i} - hence the structure of the matrices A and B in (|1.4[) . 

The discretization (|1.5I) of the linear boundary condition (ll.2j) might be 
interpreted as a Dirichlet-type condition, since the relevant subsystem of ODEs 
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involving the matrix C is easily solved exactly (cf. (|2.3I) below). We emphasize, 
however, that the objective of this paper concerns comparing the numerical 
solution to the exact solution of the Black-Scholes PDE with the linear boundary 
condition (|1.2[) at the upper boundary s = S, and not with a Dirichlet condition. 
As it turns out, the analysis in the present situation, of (jl.2p . encounters a 
variety of additional difficulties. 

Our analysis commences with an investigation of the stability of the FD 
discretization (ll.3[) . (|1.4p . This pertains to the derivation of rigorous bounds on 
the norm of the matrix exponential e tM . We shall deal here with the maximum 
norm. By | • |oo and || • ||oo we denote the maximum norm of real vectors and 
matrices, respectively. An important tool is the logarithmic maximum norm, 
which is defined for any square matrix X by 

Wi + txw^-i 

/•too A = hm , 

tj.o t 

where / is the identity matrix of the same size as X. Upon writing X = (£y )7\-=i 
a convenient formula for the logarithmic maximum norm is 

fJ>oo[X] = max {ta + V fa \ }. (1.6) 

l<z<m * — ' 

A key property is given by the following theorem; see e.g. [3) IT} l9| ITT] . 

Theorem 1.1 Let wel. Then: ^[X] < lo ||e tx ||oo < e tw (/ > 0). 

We note that we previously used the logarithmic norm in analyzing the stability 
of discretizations of the Black-Scholes and Heston PDEs when provided with 
Dirichlet boundary conditions, sec 5, 6, 12 . 

An outline of the rest of the paper is as follows. 

In Section [2] we investigate the stability of general semidiscretizations (|1.3j) . 
(|1.4|) of the Black-Scholes PDE with the linear boundary condition. We prove 
sharp upper and lower bounds for 1 1 e*^ x 1 1 oo ■ 

In Section [3] various well-known FD discretizations are considered. For each 
discretization a practical sufficient condition is obtained such that the stability 
result of Section [2] holds. 

In Section |4] we derive a convergence estimate for general semidiscretizations 
(|1.3[) . (|1.4|) of the Black-Scholes PDE with the linear boundary condition. In 
Section [5] it was found that ||e tM ||oo is essentially inversely proportional to the 
mesh width h m+ 2- We prove however the positive result that this growth, as 
h m +2 tends to zero, generally has no adverse effect on the convergence behavior. 

In Section [5] extensive numerical experiments are presented regarding the 
stability and convergence results of Sections [21 01 

In Section [5] we consider the discretization in time and prove stability and 
convergence results for the popular family of ^-methods. These results can be 
regarded as analogues of those obtained for the semidiscretization. 

In Section [7] conclusions and issues for future research are given. 



4 



2 A general stability theorem 



In this section we consider general matrices M of the form (|1.4p and derive a 
useful inclusion for the maximum norm of e tM for t > 0. We start with three 
lemmas. 



Lemma 2.1 If ZioZds i/iai 

AM _ ( etA I r 

[o 

Proof Consider the system of ODEs 

U'(t) = MU{t) 

with solution given by 

U(t) = e tM U(0). (2.1) 
Let the vector U(t) be splitted into two parts, 

' V(t) ' 



m ~ V w(t) J > 

where V(t) is an m-vector and W(t) is a 2-vector. In view of (jl.4l) one has 

fV(t) = AV(t)+BW(t), 
= CW(i). 

Thus W(t) = e* c VF(0) and 

T/'(r) - AV(t) = Be TC W(0) 

A (e- rA F(r)) = e- rA Be rC ^(0) 
ar 

e^V^t) - 7(0) = / e- TA Be TC W(0) dr 
Jo 

V(t) = e tjl V(0) + / e(*- T ^Be TC dr W(0). 
./o 

Comparing with (|2.1I) . the result of the lemma is obtained. 



□ 



The next lemma gives the maximum norm of e tC . 
Lemma 2.2 It holds that 

\\e tC \\ 00 =e-r t + (l-e-*)-^-. (2.2) 
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Proof The two eigenvalues of C are and — r with corresponding eigenvectors 
(s rn+1 S) T and (1 1) T . Thus 



which gives 



s m +i I \ 1 
S 1 \ e"' 



■Sm + l 1 

5 1 



Sm+l s m+x (l - e rt ) 



(2.3) 



Hence, 



„tC| 



max{|^e-'' t - s m+1 \ + s m+1 (l - e - rt ) , S(l - e~ rt ) + S- s m +ie- rt } 



It is readily seen that 

\Se~ rt - s m+ i\ + s m+ i(l - e- rt ) < S{1 - e- rt ) + S- s m+1 e 
Therefore, 



5(1 - e- rt ) + S- s m+1 e 



L m+2 



e- rt + (1 - e- rt ) 



2S 



hn+2 



□ 



Lemma 12.21 shows that for any given t,r,S > the maximum norm of e tC 
is essentially inversely proportional to the mesh width h m+ 2- The growth of 
He^Hoo as h m+2 decreases corresponds to the fact that at the grid point s = S 
a downwind scheme is used for the advection term in the Black-Scholes PDE. 

Let e m denote the m-dimensional unit vector (0, . . . , 0, 1) T . 

Lemma 2.3 If A is invertible, /j.<x>[A] < and a m + \/3 m \ < 0, then 

1 



| A 6 m | oq 



C^m | @m | 



Proof (i) Consider first /^[A] < 0. Define v = A 1 e m , i.e., Av — e m . Writing 
v = (vi, V2, ■ ■ ■ , v m ) T this gives the system of equations 

aivi + 7i« 2 = 0, 

favi-i + aiVi + 7ifi+i =0 (2 < i < m — 1), 

We prove by induction that \v\\ < \v2\ < •■• < \v m \. In view of (ll.6[) . the 
assumptions on A imply that all on < and ct\ + (71 1 < 0. Using this, yields 



7i_ 



1^2 1 < |u 2 |. 



G 



Next suppose |i>i-i| < \vi\ for some 2 < i < m — 1. We show that |uj| < |t>i+j.| 
By (jl.6p there holds |/3i| + CKj + |7i| < and consequently 



> and -, — - ' - ; < 1. 



We have 



+ ajUi + 7i«i+i = 
=> I a* I • |«»| < | A | • + |t* I • 
=> K| ■ \Vi\ < \/3i\ ■ \vi\ + |7i| • 

Kl - \Pi\ 

\Vi\< \Vi+l\. 

This proves the induction step, and it follows that |A _1 e m | 00 = |«| 
Subsequently, 

^mVm 1 Pm^m—l 
=>■ | "ml ' \Vm\ < 1 + 1/3 m * — 1 1 

l 



|a m | - |/3m| 

where it is used that \a m \ — \Pm\ > 0. Since \a m \ = —oe m the bound of the 
lemma is obtained. 

(ii) Consider next ^too[^4] = 0. Define the matrix A e = A — el with e > 0. 
There holds /^oo[-<4 E ] = — £ < and thus we can apply the result from part (i): 

\A S Cm|oo 



—a m + e- \p m \ 
Taking the limit e J, in this inequality completes the proof. 

□ 

The following theorem is the first main result of this paper. It provides a 
tight inclusion of the maximum norm of e for matrices M of the form (ll.4p 
and reveals that this norm is essentially inversely proportional to h m+ 2- 

Theorem 2.4 If 

rl + A is invertible, Hoo[rI + A] < 0, r + a m + \(3 m \ + \j m \ < 0, (2.4) 
then 

e- rt + (1 - e" rt ) < We^Woo < + (l + Ze^) . (2.5) 

h-m+2 h-m+2 
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Proof We employ the formula for e tM given by Lemma \2. II First, it is clear 
that | |e* A ' J | |oo > | |e tc, | |oo and by Lemma l2~2"l the lower bound is directly obtained. 
To prove the stated upper bound, we consider the maximum norm of 

' e^ A Be TC dr. 



Using formula (|2.3[) . it is readily shown that 

Be T ° = — — — [ (Se^ rT - s m+ i)e m {s m+1 - s m+1 e^ rT )e m ] . (2.6) 

For any given real numbers (j)o, 4>\ consider the vector 

/(r) = (0 o + ie ^ r )e m . 
A straightforward computation yields 



i 

o 



^ A f(r)dr = M? tA ~ I)A- l e m + ^{e tA - e - rt I)(rI + A)~h 



Note that [ioc>[rI + A] < means /ioo[A] < — r. In view of this and the assump- 
tions of the theorem it holds that both A and rl + A are invertible and, by 
Theorem O 

||e tA |U <e~ rt . (2.7) 
Consequently, we have the bound 

I f ef-^/Mfr^ < (H- e ^ t )-|^o|-|^ 1 e m |oo + 2e- rt -| ( /» 1 |-|(r/+A)- 1 e m | 00 . 
Jo 

Now let /(t) represent any column vector of Be rC . It is clear from (|2.6[) that 
|0o|<l7m|-T an d \(fri\ < |t« ' 



hn+2 n m +2 

If 7m = 0, then the result of the theorem is obvious; in fact ||e ||oo = ||e tc, ||oo 
in this case. Thus assume j m ^ 0. Then r + a m + \(3 m \ < and application of 
Lemma 12.31 to both A and rl + A yields 

\A~ 1 e m \ oc < _ and \(rl + A)~ 1 e m \ 00 < — — — 1 — — 

<^m | Pm | T Qtfn | Pm \ 

Since 

l7m| <1 and l7m| - , < 1. 



-a m - |/? m | " -r - a m - |/3 m | 

it follows that 

1 e^ A f(r)dr U < (1 + ier rt ) ■ 

"m+2 
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Using that Be rC has two columns, we arrive at the bound 

|| /* e^ A Be rC dr |U < (1 + 3e~ rt ) • . (2.8) 

JO "m+2 

Finally, in view of Lemma 12. 1\ 

~' M " <max/||e M l| 00 +ll / e^ A Be rC drV 



1 1 e | |oo s max< ||e ||oo + || y e v ' ±fe dr ||oo , ||e 

By invoking (|2.7|) , (|2.8|) and Lemma 12.21 the upper bound of the theorem is 
obtained. 

□ 

In Section[31 applications of Theorem l2.4l to various actual FD discretizations 
of the Black-Scholes PDE with the linear boundary condition shall be discussed. 
In Section@]the stability results from the present section shall effectively be used 
in the convergence analysis of FD discretizations. 

Remark 2.5 A more direct way to arrive at a bound for ||e' M ||oo is by using 
the inequality 

e (t-T)A Be rC dT n < /"|| e (*-r)X|| . M^, . || e rC|| rfr 



and then applying &2.fy) , \2. 7| ) and calculating the obtained (simple) integral. 
However, this leads to an upper bound which is substantially less favorable than 
that of Theorem \2.4\ The reason for this lies in the fact that it yields a factor 
\lm\/h m +2 and in oU actual applications the quantity \j m \ is itself inversely 
proportional to one or more mesh widths, cf. Section^ 

In subsequent applications of Theorem 12.41 the following lemma is useful. 
Lemma 2.6 If A satisfies the conditions 

' > (2 < j < m), 

7i>0 (1<J<to-1), 
ai +7i < 0, 

aj+Pj+lj = (2<j<to-1), 
a m + (3 m < 0, 

then A is invertible and fioo[A] = 0. 

Proof The conditions of the lemma directly imply, by (|1.6p . that /ioo[A] = 0. 
We next show that A is invertible. Upon setting f3\ = —u\ — 71 > and 
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7m = — OL m — j3 m > and using that etj + ft + 7 j = wc obtain 



/ ai 71 

ft Q?2 72 



\ 



Pm— 1 &rn— 1 



V 

/ -ft 71 

-ft 72 



/ 1 



m 



■Pm-1 7m- 1 



1 1 



\1 1 



=: P. 



1 7 



-7n 



Next, modify the matrix P by subtracting column 1 from columns 2,3,..., m— 1. 
This leads to 



/ -ft -ai ft /? 




Clearly, if P is invertible, then so is A. We prove that 

x e M m , Px = =>> x = 0. 

Write x = (xi, X2, . . . , a; m ) T . Then = is equivalent to the system of 
equations 



-ftxi — aia; 2 + ftx 3 + ftx 4 + . . . + ftx m _i 

- ftx 2 + 72^3 

- ^3^3 + 73^4 



= 
= 
= 



(2.9) 



-7mXl 



Pm—X^m—\ "T" 7m-l^m 



We distinguish three cases. 

(a) Assume ft ^ whenever 1 < j < m. Then ix, < 0, ft > 0, 7j > for 
all j. Starting with the last equation of (|2.9I) and moving upwards, one finds 
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that 



£0 ~\_ ■ i 



n 

k=j 



Pk J Ct r , 



■ X\ 



(2<j<m - 1). 



Substituting this into the first equation of (12.91) yields 



-j8i - ai I] 



,fe=2 



/3/c / a„ 



m—1 I m — 1 

■AE n 

J=3 \ fc=J 



7fc I 7f> 



xi = 0. 



It is easily seen that the coefficient of x\ in the latter equation is nonzero. Thus 
Xx = 0, and consequently x = 0. 

(b) Assume /3i = 0. Since ai and all jj are nonzero, (|2.9|) directly yields 
that x = 0. 

(c) Assume /3j = for certain 2 < j < m. This induces a natural partitioning 
of the matrix A where each diagonal block belongs to either case (a) or case (b) 
above. Using this, it readily follows that if Ax — then x = 0. 



□ 



Remark 2.7 Under the assumptions of Lemma \2.6\ it holds that —A is a so- 
called M-matrix. This follows, e.g., by using condition (M37) in JU Chapter 6J. 

3 Application of the general stability analysis 

In the following we apply the general stability analysis of Section [5] to actual FD 
discretizations of the Black-Scholes PDE with the linear boundary condition. 
Pertinent to the advection term on [s\, s m ] we consider three FD schemes: 

w u(s J+1 ,y u ( S] ,t) ^ (3 ia) 

hj+i 

u s ( Sj ,t) « »(^ + i,*)-»(^-i.*) > (3 . lb) 
~ -^ufa-i,t) + hi +]~ Hj u( Sj ,t) + —^L- u (s j+1 ,t). (3.1c) 

The scheme ()3.1r ) has a first-order truncation error and is called the first- order 
forward scheme. The scheme p. lb ) possesses a second-order truncation error 
whenever the grid is smooth; for the scheme p. lb ) this holds for arbitrary grids. 
We refer to (|3.1b ) and p. lb ) as the central scheme A and the central scheme B, 
respectively. Notice that these two schemes are identical if the grid is uniform. 
For the diffusion term on [s 1; s m ] we use the standard central FD scheme 

2 2 2 

u M (sj,t) « — — -u(sj-i,t) - — u{sj,t) + —u(s ]+1 ,t) , (3.2) 

hjHj tijhj+i n j+1 Hj 



11 



which possesses a second-order truncation error on smooth grids. 

Based on (I3.1[) . (|3.2I) we consider five well-known FD discretizations of the 
Black-Scholes PDE. These discretizations differ in the treatment of the advec- 
tion term rsu s on the interval [si, s m \. The linear boundary condition is always 
discretized on [s m +i, Sm+2] according to (ll.5p . For each semidiscretization, suc- 
cessive application of Lemma |2"1)1 (with A replaced by rl + A) and Theorem l2.4l 
directly gives a condition on r, a and the grid such that the stability bound 
(1231) holds. 

1. Forward. Using the first-order forward scheme give^l for 1 < j < m, 

,2 



7j = 



,2 



hjh 



jUj+l 



lj+l 



a 



lj+l 



(3.3) 



The bound (|2.5p holds for all r > 0, a > and all grids. 
2. Central A. Using the central scheme A gives for 1 < j < m, 



ft 



H = 



+ a 2 ^- 

h j Hj 

2 

2 s j 
r — — ■ 

hjhj+i 
hj+iHj 



(3.4) 



The bound holds if for all 1 < j < m: 



< r < -^-cr 2 . 



(3.5) 



3. Central B. Using the central scheme B gives for 1 < j < m, 

p . ,3 ^+1 



(7 



S 2 
2_Zj_ 

hjHj 



a j 



. s j( h 



j'+i 



5 j j 
ij+iHj 



h j+1 Hj 



The bound (f2~5|) holds if for all 1 < j < m: 



hjh 



(3.6) 



< r < t-^—o- 2 . 

hj+i 



1 Note that 0x need not to be defined, but it is convenient for the analysis. 



(3.7) 
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4. Mixed A. This scheme is defined as a suitable combination of Forward 
and Central A above. For each given 1 < j < m: if (|3 . 5[) is fulfilled, then 
define ctj, (3j, -fj by (|3.4j) : else by (13.31) . This scheme has been considered 
e.g. in [T3]. Clearly, the idea is to switch from the central scheme (|3.4p 
to the forward scheme (|3.3p in those grid points that would give rise to a 
strictly positive logarithmic maximum norm of rl + A. By construction, 
it follows that the bound (|2.5[) holds for all r > 0, a > and all grids. 

5. Mixed B. This scheme is defined as a suitable combination of Forward 
and Central B above. For each given 1 < j < m: if (|3.7[) is fulfilled, then 
define ctj, f3j, jj by (|3.6p ; else by (|3.3I) . Again, it follows that the bound 
(|2.5p holds for all r > 0, a > and all grids. 



4 A general convergence result 

In this section we prove a convergence result for general semidiscretizations (|1.3[) . 
p.4p of the Black-Scholes PDE with the linear boundary condition. Let U be 
the exact solution to (|1.3p . (|1.4[) and, for < t < T, let be the vector of 

the same size as U (t) given by 

Uh(t) = (u(si,t),u(s2,t), . . . ,u(s m +2 1 t)) T , 

where u is the exact solution to the initial-boundary value problem for the 
Black-Scholes PDE (|l.ip on < s < S with linear boundary condition (|1.2|) . 
Define the spatial discretization error 

e h {t) = u h {t)-U(t) 

and the spatial truncation error 

5 h {t) = u' h (t) - Mu h {t) - b(t). 

A standard approach to estimate spatial discretization errors is to combine an 
estimate for the spatial truncation errors with a stability bound, cf. e.g. [7]. 
However, a direct use of the bound on ||e tM ||oo from Theorem 12 .41 does not lead 
to an optimal result. To obtain a useful result in the present situation where the 
linear boundary condition is employed, we consider a partitioning of the spatial 
truncation error vector into two parts, corresponding to the intervals [si,s m ] 
and [s m+ i,s m+2 ]: 

Sh(t) = ( fx® j with 5k{t) e R m and (t) 6 R 2 . 

Using the individual stability bounds derived in Section [21 we have as a prelim- 
inary result 
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Lemma 4.1 Assume (|2.4I) holds. Then the spatial discretization error satisfies 
|e*(<)|oo<t- max {|#(0)|oo + t^- • I^Wk) (t>0). 



Proof From 
one has 

and, since £^(0) = 0, 



u' h (t) = Mu h (t) + b(t) + S h (t), 
U'(t) = MU(t) + b(t) 

e' h {t) = Me h {t) + 5 h {t) 



Lemma 12.11 yields 



£h(t) = 






Jo 




ft-» e (t-#-T)A Be TC dT 


o 





By (|2.7p . (I2~8t and Lemma [2T2| the following bounds hold whenever < <d < t: 

\\e {t ~ €)A \\oo < 1, 



Jo 



h m +2 
2S 



h 



m+2 



Hence 



Together with the integral representation above, this gives the bound on the 
maximum norm of eh(t). 

□ 

The following theorem forms the second main result of this paper. It gives 
a useful convergence estimate for general semidiscretizations (|1.3|) , (|1.4j) of the 
Black-Scholes PDE with linear boundary condition. Its proof is obtained by 
combining Lemma |4 . 1 1 with a bound for 8^. 

Theorem 4.2 Let h* > be given and assume that on [S — h* , S] x [0, T] the 
partial derivative u SS s exists and is continuous. Define 

K = 4a 2 S 3 + 6rS 2 h* and r)(t) = max{|u sss (£, t)\ : S - h* < £ < S}. (4.1) 

Assume (|2.4p holds. Then the spatial discretization error satisfies 

\sh{t)\oo < i-max {l^^loo + nnii})} whenever < h m+2 < h* , < t < T. 
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Proof Write s = s m+ i and h = h m+ 2- Pertinent to the point s m+ \ we have 

u(S, t) — u(s, t) 



5*(t) = u t (s,t) 



h 



* l a 2 s 2 u ss (s,t) + rs 



u s (s,t) 



+ ru(s, t) 

u(S, t) — u{s, t) 



By Taylor's theorem and using u ss (S,t) = there follows 

u(s,t) = u(S,t)-hu s (S,t) + E (h,t), 
u s (s,t) = u s (S,t)+E 1 (h,t), 
u ss (s,t) = E 2 (h,t), 

where 

E (h,t) = -lh 3 u sss (£ ,t) , = \h 2 ,t) , E 2 (h,t) = -hu sss (&,t) 

with certain £o,<!;i,£2 in (s,S). Substitution into the above formula yields 
S^it) = ±a 2 s 2 E 2 (h,t)+rsE 1 (h,t) + ^-E (h,t), 



h 



which readily leads to the estimate 



K 1 {t)\<{¥ 2 S 2 + lrSh}i 1 (t)-h. 



Analogously, pertinent to the point s m +2, there holds 



/c ,x c u(S,t) - u(s,t) 
Ut[o, t) — rb h ru(b, t) 



rS 



u B (S,t) - 



u(S, t) — u{s, t) 



r 9 

^-E (h,t) 



and 

It thus follows that 
and 



\^. 2 (t)\ < k rSr l(t) ■ h 2 



\^{t)U<{ 2 -a 2 S 2 + lrSh*} V {t)-h 



\ ■ |^(t)|oo < «T?(*). 



(4.2) 



Application of Lemma |4~T1 then gives the desired estimate for \eh{t)\oo- 



□ 
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The estimate of Theorem 14.21 for the spatial discretization error consists of 
two contributions, corresponding to the two intervals [s\, s m ] and [s m +i, s m+ 2\- 
The first contribution is equal to the part of the spatial truncation error per- 
tinent to [si,s m ]. For any given FD discretization this can be estimated in a 
standard way by Taylor expansion. The second contribution is equal to Kr)(-&) 
and depends on the partial derivative u sss of the exact option price near s = S. 
For a wide range of financial options it is plausible that this contribution can 
be made arbitrarily small upon taking the upper bound S sufficiently large (in 
the case of European call and put options this is readily proved). Theorem 14.21 
thus expresses the useful result that the contribution to spatial discretization 
error of the semidiscretized linear boundary condition is negligible, provided S 
is sufficiently large. Key to the proof is that the stability bounds derived in 
Section [2] admit a growth of errors from the interval [sm+i> s m+2] that is at 
most inversely proportional to h m +2 but this growth is precisely offset by the 
factor h m +2 that arises in the part of the spatial truncation error pertinent to 
this interval. 

We note that in Theorem 14. 21 it is tacitly assumed that the initial function 
is smooth, which has been used for ease of the analysis, but is often not fulfilled 
in applications. The numerical experiments in the subsequent section deal with 
a nonsmooth initial function. 



5 Numerical experiments 

In this section the stability and convergence results of Sections [3J 0] are illus- 
trated by numerical experiments. We consider the five FD discretizations of the 
Black-Scholes PDE formulated in Section [3] with the linear boundary condition 
discretized as given in Section [1] For each FD discretization, we also consider 
its analogue where the definition of otj , /3j , 7j for 1 < j < m in Section [3] is ex- 
tended to j = m + 1. This corresponds to the approach by Windcliff , Forsyth & 
Vetzal |J3] where in the penultimate grid point s m +\ no modified discretization 
is applied related to the linear boundary condition. In the following we shall 
refer by LBC1 to the numerical treatment of the linear boundary condition as 
defined in Section [1] and by LBC2 to the treatment as considered in [T5]. 

For the numerical experiments a typical smooth, non-uniform spatial grid is 
chosen. Let E £ (0, S) and c > be given and fixed. Consider the continuous, 
strictly increasing function 

tp(£) = E + c- sinh(C) (a < £ < b) 

with a = sinh _1 (-,E/c) and b = sinh _1 ((S" - E)/c). Let 

£j=a + j-Af (0<7<m + 2) with Af = & - - . 

Then a non-uniform grid = So < s\ < . . . < s m+ i < s m+ 2 = S is defined by 
the transformation 

(0<i<m+2). 
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The parameter E can be viewed as the exercise price of a vanilla option and c 
determines the fraction of grid points Sj that lie in the neighborhood of E, 

hj » c • A£ whenever Sj « E. 

The grid is smooth in the sense that there exist real constants cq, c\, C2 > 
(independent of j and m) such that the mesh widths hj satisfy 

c • A£ < hj < ci ■ A£ and - hj\ < c 2 ■ (A£) 2 . 

The above type of grid is often used in financial applications, cf. e.g. [10]. We 
(arbitrarily) set E = 100, c = £/5, 5" = 400. 

5.1 Stability experiments 

First, the stability of the ten FD discretizations discussed above is numerically 
investigated. For each FD discretization we consider the maximum of ||e* M || 00 
over t > for the dimensions m = 50, 100, . . . , 1000. We employed the Matlab 
function expm to evaluate the matrix exponential and computed the values of 
\\e tM | |oo for t = 0, 1, . . . , 100; this was found to give a reliable estimat^] for the 
maximum over all t > 0. For the experiments three pairs r, a have been chosen. 
Figure Q] displays the obtained results. On the top row r = 0.1 and a = 0.3; on 
the middle row r = 0.3 and a — 0.1; on the bottom row r = 0.2 and a = 0. The 
left column represents the LBC1 treatment of the linear boundary condition 
and the right column represents the LBC2 treatment. Note that the scales on 
the vertical axes vary. 

If r = 0.1 and a = 0.3, then the Mixed A and B discretizations are identical 
to Central A and B, respectively; in this case the conditions (I3.5[) . (|3.7I) always 
hold. On the other hand, if r — 0.2 and a = 0, then Mixed A and B both reduce 
to the Forward discretization; in this case (|3.5I) . (13.71) never hold. Finally, if 
r = 0.3 and a = 0.1, then Mixed A and B form an actual mix of Forward with 
Central A and B, respectively. 

Considering the left column of Figured! a main observation is that the results 
for the Forward and the Mixed A and B discretizations with the LBC1 treatment 
agree with the stability bound (|2.5p given by Theorem 12.41 Indeed, for all r, a 
pairs the maximum of He'^Hoo over t is found to be directly proportional to m, 
and moreover, 

max He*"! |oo « . (5.1) 

t>o h m+2 

If r = 0.1 and a = 0.3, then the results for Central A and B with LBC1 clearly 
also agree with (|2.5p . as these discretizations coincide with Mixed A and B. We 
do not have a theoretical stability result for the central discretizations for the 
other two r, a pairs, but it is interesting that the estimate (|5 . 1 [) is also observed 
in the left column for Central A and B if r = 0.3 and a = 0.1, and for Central A 

2 Except for Central A and B with LBC2 if <r = 0, see the discussion on this case in the 
text further on. 
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LBC1, r = 0.1,0- = 0.3 



LBC2, r = 0.1, a = 0.3 



1000 




1000 



6.5 



5.5. 
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-e- Mixed B 
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LBC2. r = 0.3, o — 0.1 



1000 



120 




100 
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LBC2. r = 0.2, a = 
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i 
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-B- Mixed A 
-e- Mixed B 
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Figure 1: Plot of max{||e 



tMl 



:t = 0,1,..., 100} vs. m = 50, 100, 



, 1000. 

On the top row r = 0.1 and a = 0.3; on the middle row r = 0.3 and a = 0.1; 
on the bottom row r = 0.2 and a = 0. The left column concerns the LBC1 
discretization of the linear boundary condition; the right column concerns the 
LBC2 discretization. In all cases E — 100 and S — 400. Notice the different 
scales on the vertical axes. 
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if r = 0.2 and a — 0. Although (15.11) is not found for Central B in the latter case, 
the maximum of He^Hoo over t still appears to be at most directly proportional 
to to. 

Considering the right column of Figure [TJ we observe in the top and middle 
rows that for all five FD discretizations with the LBC2 treatment there appears 
to be an upper bound on ||e tM ||oo that is uniform both in t and the dimension to. 
This upper bound is small in the top row and larger in the middle row. It appears 
to increase if the ratio a 2 /r decreases. In the limit case of a = 0, displayed in 
the bottom row, one obtains again linear growth with to for the Forward and 
Mixed A and B discretizations with LBC2. In fact, it is readily seen that these 
three discretizations then all reduce to Forward with LBC1. We notice that 
similar observations as here concerning LBC2 were made in [14]. 

If a = 0, then the results for the Central A and B discretizations with LBC2 
are not displayed in the figure. For each considered dimension m, the obtained 
maximum of He^Hoo is very large in this case. For example, if m = 100, then 
the maximum values (over t = 0, 1, ... , 100) for Central A and B are equal to 
2.5 x 10 3 and 1.5 x 10 5 , respectively, and if to = 200, then they are equal to 
9.8 x 10 3 and 5.8 x 10 5 . In fact, there appears to be growth with to 2 . 

When a is not small, then comparing the results in the left and right columns 
of Figure Q] one may be inclined to prefer the LBC2 treatment over LBC1: for 
both r, cr pairs with nonzero a one observes, for each given FD discretization, 
that the maximum of ||e tM ||oo is more favorable for LBC2 compared to LBC1. 
However, considering the actual convergence behavior of the FD discretizations, 
we find no essential difference between the two treatments, as is illustrated next. 

5.2 Convergence experiments 

Here we numerically examine the convergence behavior of the ten FD discretiza- 
tions of the Black-Scholes PDE with linear boundary condition discussed above. 
For the experiments we consider a European call option with exercise price E, 
so that 

u(s, 0) = max(0, s — E) (s > 0), u(0,*)=0 (0<t<T), 

and set E = 100, T — 5, S = 2000. In the experiments we compute the maxi- 
mum norm of the spatial discretization errors 

e h (T) = u h (T)-U(T) 

for a sequence of values to with 10 2 < to < 10 4 . Here Uh denotes the restric- 
tion to the spatial grid of the exact call option price function, given by the 
Black-Scholes formula. The semidiscrete solution vector U(T) to (1 1 .3[) . (|1.4I) is 
approximated, with sufficiently high accuracy, by applying the Crank-Nicolson 
method using N = 10 4 time steps. In view of the nonsmooth initial condition, 
two initial (damping) substeps are taken with the implicit Euler method. This 
approach is often referred to as Rannacher time stepping. 
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LBC1, r = 0.1,<7 = 0.3 



LBC2, r = CI, <r — 0.3 




Figure 2: Plot of \eh(T)\oo vs. 1/m for values 10 2 < m < 10 4 equally spaced 
on a logarithmic scale. On the top row r = 0.1 and a = 0.3; on the bottom row 
r = 0.3 and a = 0.1. The left column concerns the LBC1 discretization of the 
linear boundary condition; the right column concerns the LBC2 discretization. 
In all cases E = 100, T = 5, S = 2000. 
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The results are displayed in Figure O On the first row r — 0.1, a = 0.3 
and on the second row r — 0.3, a = 0.1. The left column represents the LBC1 
treatment of the linear boundary condition and the right column represents the 
LBC2 treatment. 

As a first observation it is clear that LBC1 and LBC2 almost always lead to 
the same spatial discretization error in the experiments. 

If r = 0.1 and a = 0.3, then the Mixed A and B discretizations are identical 
to Central A and B, respectively. For these four discretizations we find, by least- 
squares approximation in the region to < 5000, an order of convergence equal 
to 2.0. Once to gets larger than 5000, then the (fixed) error due to the linear 
boundary condition dominates. The Forward discretization has an observed 
order of convergence equal to 1.0. 

If r — 0.3 and a = 0.1 the observed orders of convergence are 0.9 for Forward, 
2.0 for Central A and 1.9 for Central B. The Mixed A and B discretizations are 
in this case an actual mix of Forward with Central A and B, respectively. Ta- 
ble Q] gives for each to under consideration the fraction of grid points where the 
Forward discretization is used. Clearly, for small m this fraction is large and 
Mixed A and B show spatial discretization errors similar to those for Forward, 
whereas for large to this fraction is small and Mixed A and B show errors simi- 
lar to those for the Central A and B discretizations. A rigorous analysis of this 
correlation will be left for future research. 



TO 


Fraction Mixed A 


Fraction Mixed B 


100 


57.0% 


58.0% 


129 


47.3% 


48.8% 


167 


34.1% 


36.5% 


215 


9.8% 


11.2% 


278 


7.9% 


7.6% 


359 


6.4% 


6.4% 


464 


5.2% 


5.2% 


599 


4.2% 


4.2% 


774 


3.4% 


3.4% 


1000 


2.7% 


2.7% 


1292 


2.1% 


2.1% 


1668 


1.7% 


1.7% 


2154 


1.3% 


1.3% 


2783 


1.0% 


1.0% 


3594 


0.8% 


0.8% 


4642 


0.6% 


0.6% 


5995 


0.5% 


0.5% 


7743 


0.4% 


0.4% 


10000 


0.3% 


0.3% 



Table 1: Fraction of grid points where the Forward discretization is used in 
Mixed A and Mixed B vs. to in the case of r — 0.3, a — 0.1 and LBC1. 
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6 Time discretization 



In this section we study the time discretization of the semidiscrete system (|1.3p . 
(|1.4I) by the well-known family of ^-methods, which includes the popular Crank- 
Nicolson method (trapezoidal rule) and implicit Euler method as special cases. 
As noted in Section [TJ the subsystem of ODEs involving the matrix C could be 
solved exactly, but it is more interesting and useful, both from a theoretical and 
practical point of view, to consider the time discretization of the semidiscrete 
system (|1.3p . (|1.4[) as a whole. 

Let parameter 9 s [\, 1] be given and fixed. Let step size At = T/N with 
integer N > 1 be given and define t n = n At, b n = b(t n ) for n = 0, 1, 2, . . . , N. 
The 6 '-method generates, in a successive way, for n = 1, 2, . . . , N an approxima- 
tion U n to U(t n ) by 

U n = U n -x + (1-0) At [MU n -i + bn-!) + 9 At (MU n + b n ). 

The choices 9 — | and 9 = 1 yield, respectively, the Crank-Nicolson method 
and implicit Euler method. The above recurrence relation can be written as 

U n = (p(AtM)U n -i + (J - 9 At Af)" 1 [(l - 9)Atb n - 1 +9Atb n ], (6.1) 

where <p is the so-called stability function of the method, given by 

l + (l_0)z (forz£C)! 
1 — 9z 

and 

(p(X) = (I - 9X)-\I + (1 - 9)X) = (/+(!- 9)X)(I - 9X)- 1 



for square matrices X such that I — 9X is invertible. 

We first study the stability of the fully discrete process (|6.1I) . The subsequent 
two lemmas can be viewed as analogues of Lemmas l2.1[ 12.21 for the semidiscrete 
system. 

Lemma 6.1 For 1 < n < N there holds 
<p(AtM) n 



<p(At A) n 




O 


cp(AtC) n 



wher 



Yi 



(I - 9AtAy 1 AtB{I -9AtC)-\ 



Y n = Y,<p(AtA) n -i- 1 Y 1 <p(AtC) j . (e 
Proof The formula is readily obtained by induction to n and noting that 



(6.2a) 
(6.2b) 



(I -9 At A)- 1 


9Y 1 


O 


(I- 


-9AtC)~ L 



(6.3) 
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□ 



Lemma 6.2 For 1 < n < N there holds 



MAtcy 



n 



x n + (l - x n ) 



2S 



with x — ip(~rAt). 



h 



<m+2 



Proof With the eigendecomposition of C it is easily verified that 



(f(AtC) 



Tl 




Sx n - s m+x s m+i (l-x") 
S(x n - 1) S - s m+1 x n 



) 



and the rest of the proof is similar to that of Lemma l2~2l using \x\ < 1. 



Concerning the discretization on the spatial domain [s%, s m ] we shall assume 
in the following that there exists a real constant K , independent of the dimension 
m and number of time steps n, such that 



\\<p(At A^Hoo < K whenever At = T/N, 0<n<N, N>1. (6.4) 



In the literature much attention has been paid to establishing (|6.4[) , under a va- 
riety of conditions on the matrix A. For the implicit Euler method (6 = 1) the 
neat result is well-known that (|6.4|) is fulfilled with K = 1 whenever /ioo[^] < 0, 
cf. e.g. [4J [8] . Hence, this is guaranteed under the condition (|2.4[) . For all other 
time discretization methods, however, the available results in the literature im- 
plying (16.41) with a constant K independent of m and n require, to the best 
of our knowledge, stronger conditions on A. Notably, conditions on the resol- 
vent, the numerical range and pseudospectra have been extensively investigated 
in the literature, cf. e.g. [9j [11]. Some of these results appear to be useful in 
our current application, but a verification of the pertinent conditions on A is 
highly non-trivial. As our main interest in this paper lies in studying (the im- 
plications of) the discretized linear boundary condition on [s m+ i, s m+ 2\, which 
corresponds to the matrices B and C, we shall leave the analysis of (|6.4[) when 
i < 9 < 1 for future research. 

The next theorem can be regarded as a discrete analogue to Theorem 12.41 
Theorem 6.3 If {220) and (El) then forl<n<N, 



with x = (p(—rAt). 

Proof For any integer j > 0, let the rational function ipj be defined by 



□ 



x n + (l-x n ) 
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< IMAtMHoo < K + 




tpj(z) 



<p( Z y 



(z e C). 



l-6z 
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Consider the formula for <p(AtM) n given by Lemma T6. 11 The lower bound on 
its maximum norm is clear by Lemma 16.21 To prove the upper bound, write 

n-l 

Y n = <p{&tA) n -i- 1 {I-dMA)- 1 MB{I-dMC)- 1 tp{MCy 

3=0 
n-l 

= v(AtA) n - j - 1 (I - OAtA^AtBilijiAtC). 

3=0 

It holds that 

Bipj(AtC) = 7m [ (Sxj - s m+1 )e m (s m+1 - s m+1 Xj)e m ] (6.5) 

"m+2 

where Xj — ipj(-rAt) £ [—1, 1]. Both columns of this matrix are of the form 

fj = (00 + 4>\Xj)e m 

with real numbers <j>o, (f>\ independent of j. 
We have 

n-l 

At fi^tA) 71 ^- 1 (I - 6 At A)- 1 

3=0 
n-l 

= At v{AtA) k (I - OAtAy 1 

= At(<p(AtA) n - I)(<p(AtA) - I)' 1 {I -9 At A)~ l 

= At (ip{At A) n - I){I + (1 - 9) At A- I + 9 At A)- 1 

= (<p(AtA) n - I)A~ l . 
By similar algebraic manipulations, there follows 

n-l 

At ^{AtA)"-^ 1 {I - 9AtA)- 1 x J 

3=0 

= (<p(AtA) n - v (-rAt) n I){rI + A)' 1 . 

Consequently, 

n-l 

At J2 f{At A)"^ 3 ^ 1 (I - 9AtA)- 1 f j (6.6) 

3=0 

= 0o fa(Ai A) n - I)A- X e m + 0! (<p(AtA) n - <p(-rAt) n I)(rI + Ay 1 e m . 
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Continuing from here along the same lines as in the proof of Theorem 12.41 and 
using the condition (|6.4[) . we arrive at 

4(K + 1)S 

1 1 1 n 1 1 oo _ t 

n m +2 

Together with Lemma [6.2[ the stated upper bound on ||y>(At M) n \ |oo now di- 
rectly follows. 

□ 

In the subsequent convergence analysis of the process (|6.1[) . the matrix 
tp n (AtM) = <p(AtM) n (I -OAtM)- 1 
arises. By Lemma T6. II and formula (|6.3I) . 



i/j n {AtM) = 



1p n (At A) 


6 (p(At A) n Yi + Y n (I - 6 At C)- 1 








For the analysis below we need upper bounds on the maximum norms of the 
constituent sub-matrices. Putting Yq — O, there holds 

Lemma 6.4 If and flO]) then for < n < N - 1, 

HV>n(AtA)|| 
ll^n(A«C)|| 
IMAM)" Fill 

WYnii-eAtc)-^ 

Proof The bound (|67fh ) follows directly from ([6~4"]) and the fact that QU) 
implies - 9 At A)~ 1 \\ oc < 1. The bound (|67fb ) is obtained using the same 
arguments as in the proof of Lemma 16.21 In order to prove (|6.7b l we note that 

(I - 6 At A)- 1 = (<p(AtA) ~ I) (At A) -1 

and using this gives 

ip(At A) n Y\ = {ip{AtA) n+1 -ipiAtA^A^Bil -OAtC)- 1 . (6.8) 

By (US]) with j = 0, 

B (I - OAtCy 1 = — — — [ (,5*a;o - s m+ i)e m (s m +i - s m+1 x )e m ] 

"m+2 



< K , 



< 



< 



< 



AS 

h m +2 

AKS 

h m +2 

A(K+l)S 
h m +2 



(6.7a) 
(6.7b) 
(6.7c) 
(6.7d) 



25 



where xq = (1 + OrAt)^ 1 e (0, 1). If /o represents any of the two columns 
of this matrix, then by a same argument as in the proof of Theorem 12.41 there 
follows 

L4-7oU < j^— ■ 

Consequently, 

S AKS 



\\(ip(At A) n+1 - ip(AtA) n ) A~ X B (1 - 9AtC)~ 1 \\ 00 <2K-2- 



hn+2 "m+2 



The proof of the bound (I6.7H ) is identical to that for | \Y n \ \oo given above, except 
that ipj( z ) needs to be replaced by ij)j{z)/{l — 9z). 

□ 

To prove the convergence result for the time discretization process (|6.ip , we 
also need the following result. 

Lemma 6.5 Assume (|2.4[) and (|6.4[) hold. Let g : [0, S] — > R be any given con- 
tinuously differ entiable function and 

w = ( 9{ ( Sm+1 l ) e M 2 . 

V 9( s rn+2) J 

Then there exists £ £ (s m+ i, s m +2) such that for all < n < N — 1: 

|^„(AiCHoo < \g(S)\+2S\g'(0\, (6.9a) 
\^(AtA) n Y 1 w\ oc < 2K(\g(S)\ + S\g'(0\), (6.9b) 
\Y n (I -dAtC^wU < (K + l)(\g(S)\+2S\g'(0\)- (6-9c) 

Proof Write s — s m+ i and h = h m+ 2- Let £ S (s, S*) be such that 

g(s) = g(S) - ft S '(£). 
Then the vector w can be written as 

For any rational function ^ with V'(O) = 1 there holds 
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where x = tf>(—rAt). Application of this matrix to w, in the above form, readily 
yields 

MAtOw- ( 9(S)x-g'(0(Sx-s) \ 
^ AtC ) w -{ g(S)x-g>(0(Sx-S) J' 

Observe the important fact that there is no factor 1/h present here. 

(a) The bound (|6.9a ) is obtained upon taking tp = ip n and using |x| < 1. 

(b) By formula 

|<^(AtA)™yiw| 00 < 2K ■ \A- X B (I - OAtC)- 1 w^. 

Considering ip(z) = tpo(z) = (1 — 9z)~ 1 yields 

A~ l B (I - 9 At C)- 1 w = (g(S)x - g'(Q(Sx - a)) 7 m A~ 1 e m 

with a; = (1+9 rAt)^ 1 . As in the proof of Theorem l2.4l we have |7 m ^4~ 1 e m | o < 1 
and, together with < x < 1, there follows 

\A- 1 B(I-9AtC)- 1 w\ OQ < |<?(S)| +%'(£) |, 

which completes the proof of (|6.9b ). 

(c) By formula 

n-l 

Y n (I -9AtC)- 1 w = ^{AtAy 1 - 1 - 1 (I- 9AtA)- l AtBi> j (AtC)w 

3=0 

with ipj(z) — ipj(z)/(l — 9z). Taking tp — tpj in the above general formula, we 
get 

Bipj(AtC)w = (0o + <j>\Xj )e m , 

where 

xj = ijj(-rAt) , 0o = 7m g'iQs , 0i = 7m (g(S) - g'(0 s )- 

It is convenient to set Xj = ipj(—rAt) and 0i — 0i/(l + Or At). Then 0iXy = 
4>iXj and formula (|6.(3|) directly gives 

Y n (I - 9AtCy 1 w 



At ip(AtA) n - j - 1 {I - 6 At A)' 

3=0 

0o (f(AtA) n - I)A- 1 e m + 0! (ip(AtA) n - v (-rAt) n I)(rI + A)~ l e r 
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Using that \j m A 1 e m \ oc < 1 and \j m (rl + A) l e m \oa < 1 yields 

\Y n (I-9AtC)- 1 w\ 00 < (K + l)(\g'(0\S + \g(S)-g'(OS\) 
< (K + l)(\g(S)\+2S\g'(0\), 

which proves the bound (|6.9b ). 

□ 

As in Section |4l let the vector Uh{t) be given by 

Uh(t) = (u(si,t),u(s 2 , *),..., u(s m+2 , t)) T , 

where u is the exact solution to the initial-boundary value problem for the Black- 
Scholes PDE (fTTTT) on < s < S with linear boundary condition (|1.2p . The 
following theorem provides a useful estimate for the space-time discretization 
error, defined by 

e n = u h {t n )-U n (0<n<N). 

It essentially states that the estimate for the spatial discretization error from 
Theorem 14 . 21 remains valid after time discretization up to a c- (At) p term, where 
p denotes the classical order of consistency of the 0-method and c is a constant 
independent of the spatial grid and the time step. 

Theorem 6.6 Let p=lifh<9<l and p = 2 if 9 = i. Assume that all 
partial derivatives ofu of orders < p+2 exist and are continuous on [0, S] x [0, T] . 
Let h* > be given and let n, rj be defined by (|4.ip . Assume (|2.4|l and (|6.4[) hold. 
Then there exists a real constant c (depending only on u, S, 9 and K) such that 

|£nU <t n - IK ■ max |^(i9)|oo + (K+ max + c ■ (At) p \ 

whenever < h m+2 < h* , At = T/N, 1 < n < N, N > 1. 

Proof Let the local space-time error 6 n in the n-th step of (|6.ip be defined by 

u h (t n ) = <p(AtM)u h (t„-i) + (I- 9AtM)- 1 [{\ - 9) At 6 n _i + 9 At b n ] +?„. 
Subtracting (|6.ip from this yields 

n 

e n = ^(AtM)£ n _! +?„ = ...= <p(At M) n ~ j 5 3 . (6.10) 

i=i 

With the spatial truncation error 

r r ,s / ^(*) \ 
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as defined in Section [H one can express b(t) as 

b(t)=u' h (t)-Mu h (t)-5 h (t). 

Inserting this into the definition of 6j it readily follows that 

Sj = (I -6AtM)- 1 [(l-d)At5 h (t j - 1 ) + dAtS h (t j )] + 

(I - 6M M)~ l [uh(tj) - Uh^-i) - (1 - 9) Ai<(t,_i) -GAtu'^tj)]. 

The above formula for 5j consists of two terms, corresponding to the truncation 
error in space and the truncation error in time. In view of (|6.10[) . we shall study 
tp(AtM) n - j 5j. 

Concerning the first term, the bounds for the submatrices of ip n -j(At M) = 
<p(AtM) n ~i {I - OAtM)- 1 given by Lemma RTH directly lead to 

l^n-j(AtM) [(1 ~9) MShfy-i) + 9At5 h (t 3 )} U < 

At-{K- max . [^(*0loo + 4(2 f + 1)5 « 

I. 1=3-1,3 flm+2 1=3-1 J 

Invoking the estimate (|4. 2[) for 6jf this gives 

l^n-i(AtM) [(1 - 8) At5 h {t^x) + 0AtSh(tj)] U < 

At-lK- max |<^(^)U + >? • max \r](ti)\ c 

I 1=3-1,3 1=3-1,3 

where k — (K + 

Concerning the second term, for the function g : [0,S] — > K defined by 
<?(s) = u(s, iy) — it(s, tj-i) — (1 — 9) At u t {s, tj-i) — 9 At u t (s, t 3 ) 
standard Taylor expansion shows that 

\g(s)\ <c (Aty +1 , \g\ s )\< Cl (Atr +1 (0 < a < S) 

with 

P=l, CO = flKtHoo j Ci = |||u s tt||oo (if i < 6< < 1), 
P=2, Co = ^||wttt||oo , Ci = ^HUatttHoo (if # = |)i 

where | • |oo designates the maximum norm of a real function on [0, S] x [0, T\ 
Using now Lemma 16.51 and the partitioning 

u h (t)=l Vh ^ K | with v h (t)eR m and w h (t) £ R 2 , 

V w h(t) 
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it follows that 

| ^„_j (At M ) [u h fa ) - u h (t rf _! ) - (1 - 6) At u' h (tj—i) At u' h (tj )] | x < c (At) 
with 

c=(4K + l)co + (AK + 2)S Cl . 
Combining the above bounds, we are led to 

n 

|e„U < £>(AiM) n -' ^loo 



p+i 



" r 

< ^At- if - max |^(^)| 



+ k- max |j?(t;)|oo + c(At) p 



< t n -{K- max |(5,t(^)|oo + /? • max 77(?9) + c (At) p 

1 0<i?<t„ 0<iJ<t„ 



□ 



7 Conclusions 

In this paper we have analyzed the stability and convergence of discretizations, 
both in space and time, of the Black-Scholes PDE when it is provided with 
the linear boundary condition. This condition states that the second derivative 
of the option value vanishes when the underlying asset price gets large. For 
the space discretization we considered finite difference schemes and for the time 
discretization the well-known family of ^-methods. Concerning stability, we 
derived tight inclusions for the maximum norm of e tM and ip(At M) n where 
the matrix M represents the semidiscretized Black-Scholes PDE operator, ip 
denotes the stability function of the ^-method, time t > 0, step size At > 
and integer n > 0. The obtained inclusions reveal that these norms can grow 
essentially directly proportional to the dimension of the matrix M, i.e., the 
number of spatial grid points. We subsequently proved the positive result that 
this growth has in general no adverse effect on the convergence behavior of the 
discretizations. 

In future research, in addition to the issues already mentioned previously, 
we wish to extend the analysis to the discretization of the linear boundary 
condition that was considered by Windcliff, Forsyth & Vetzal [13] ■ Also we aim 
at studying more advanced, multi-dimensional PDEs in finance, such as the 
Hcston PDE, as well as other discretizations in space and time. 
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